 function [decoded_text,y] = Rx(r,rolloff,desired_user)
%% parameters
close all;
dm=0; %debug mode flag
fs=50E6;
Ts=1/fs;
Tsym=121.63e-9; %Symbol Timing = 121.6 ns
M=fs/(1/Tsym); %Oversampling Factor
datlent=224;
%% DEBUG - input
if dm==1
   load 'medium';
end
%% DEBUG - plot r(kTs)
if dm==1
    figure(1),plotspec(r,Ts);
end
%% filtering replicas with LPF cutoff at about 30MHz
fnyq=1/(2*Ts); lcut=19.9E6; ucut=20.1E6;
b=firpm(500,[0 lcut/fnyq ucut/fnyq 1],[1 1 0 0]);
rlp=filter(b,1,r);
rlp=r;
%% DEBUG - plot filtered r -> rlp
if dm==1
    figure(2),plotspec(rlp,Ts);
end
%% pllpreprocess.m: send received signal through square and BPF  
rsq=rlp.^2;                        % square nonlinearity
lcutA=18.9E6; lcutB=19.1E6; ucutA=20.9E6; ucutB=21.1E6;  %cutoff frequencies for BPF
h=firpm(500,[0 lcutA/fnyq lcutB/fnyq ucutA/fnyq ucutB/fnyq 1],[0 0 1 1 0 0]);     % narrow BPF center frequency at 20MHz
rp=filter(h,1,rsq);                  % filter to give preprocessed r

%% DEBUG - plot preprocessed r -> rp
if dm==1
    figure(3),plotspec(rp,Ts);
end
%% dualplls.m: estimation of carrier via dual loop structure
t=0:Ts:length(rp)*Ts-Ts;                 % time vector
mu1=.02; mu2=.002;                       % algorithm stepsizes
f0=(20E6)/2;                             % assumed freq. at receiver
lent=length(r); th1=zeros(1,lent);       % initialize estimates
th2=zeros(1,lent); carest=zeros(1,lent);
for k=1:lent-1                       
 th1(k+1)=th1(k)-mu1*rp(k)*sin(4*pi*f0*t(k)+2*th1(k));          % top PLL
 th2(k+1)=th2(k)-mu2*rp(k)*sin(4*pi*f0*t(k)+2*th1(k)+2*th2(k)); % bottom PLL
 carest(k)=cos(2*pi*f0*t(k)+th1(k)+th2(k));                   % carrier estimate
end
%% DEBUG - plot carrier tracking 
if dm==1
    figure(4),
    subplot(3,1,1), plot(t,th1)                         % plot first theta
    title('output of first PLL')
    ylabel('\theta_1')
    subplot(3,1,2), plot(t,th2)                         % plot second theta
    title('output of second PLL')
    ylabel('\theta_2')
    subplot(3,1,3), plot(rp'-carest) % plot difference between estimate and preprocesssed carrier rp
    title('error between preprocessed carrier and estimated carrier')
    xlabel('time')
    ylabel('f_0 - f')
    figure(5),
    plotspec(carest,Ts);
end
%% Demodulation with Tracked Carrier
rbase=rlp.*carest';
%% DEBUG - plot baseband rpl -> rbase
if dm==1
    figure(6),
    plotspec(rbase,Ts);
end
%% Match Filtering
L=8/2;                                  % pulse width in symbol durations
mfilter=srrc(L,rolloff,M,0);                    % receive filter H sub R
v=conv(rbase,mfilter); 
%% DEBUG - match filtering results
if dm==1
    figure(7),plotspec(v,Ts);
end
%% Downsampling - timing recovery
n=round(length(v)/M);              % number of data points
m=M;                               % oversampling factor
l=L;                               % 1/2 length of pulse shape (in symbols)

% clockrecOP.m:  clock recovery maximizing output power
% find tau to optimize J(tau)=|x(kT/m+tau)|^2 
% prepare transmitted signal
% run clock recovery algorithm
tnow=l*m+1; tau=0; xs=zeros(1,n);           % initialize variables
tausave=zeros(1,n); tausave(1)=tau; i=0;
mu=0.02;                                    % algorithm stepsize
delta=2;                                    % time for derivative
while tnow<length(v)-2*l*m                    % run iteration
    i=i+1;
    xs(i)=interpsinc(v,tnow+tau,l);           % interpolated value at tnow+tau
    x_deltap=interpsinc(v,tnow+tau+delta,l);  % get value to the right
    x_deltam=interpsinc(v,tnow+tau-delta,l);  % get value to the left
    dx=x_deltap-x_deltam;                     % calculate numerical derivative  
    tau=tau+mu*dx*xs(i);                      % alg update (energy)
    tnow=tnow+m; tausave(i)=tau;              % save for plotting
end
%% DEBUG - constellation diagram
if dm==1
    % plot results
    figure(8), subplot(2,1,1), plot(xs(1:i-2),'b.')    % plot constellation diagram
    title('constellation diagram');
    ylabel('estimated symbol values')
    subplot(2,1,2), plot(tausave(1:i-2))               % plot trajectory of tau
    ylabel('offset estimates'), xlabel('iterations')
end
%% Correlator 1
framlent=784;
preamble=letters2pam('A0Oh well whatever Nevermind');
prelent=length(preamble);
xsc=xcorr(preamble,xs);                % do cross correlation
[C,ind]=max(abs(xsc));                 % location of largest correlation...
headstart=length(xs)-ind+1;            % ...gives place where header starts
headstart=mod(headstart,framlent)+framlent;

numpre=0;
g=length(xs)-headstart;
while g>=0
    g=g-framlent;
    numpre=numpre+1;
end
%% DEBUG - Correlator 1
if dm==1
    figure(9),subplot(3,1,1), stem(preamble)       % plot header
    title('Header')
    subplot(3,1,2), stem(xs)             % plot data sequence
    title('Data with embedded header')
    subplot(3,1,3), stem(abs(xsc))            % plot correlation
    title('Correlation of header with data')
end
%% Equalization and Decision
cnt=1;
mu=.01; delayspread=7;                        % stepsize and delay delta
numtaps=4*delayspread; f=zeros(numtaps,1);    % initialize equalizer at 0
timer1=0;
timer2=0;
for j=0:1:(numpre-2)
    index1=headstart+framlent*j;
    % LMS Equalizer
    while timer1<=prelent-1 
        rr=xs(index1+delayspread+timer1-1:-1:index1+delayspread+timer1-numtaps);
        e=preamble(timer1+1)- rr*f;
        f=f+mu*e*rr';
        y(cnt)=rr*f;
        timer1=timer1+1;
        cnt=cnt+1;
    end
    timer1=0;
    index2=index1+prelent;
    % DD Equalizer
    while timer2<=datlent*3-1 
       rrb=(xs(index2+delayspread+timer2-1:-1:index2+delayspread+timer2-numtaps))';
       e=quantalph(rrb'*f,[-3 -1 1 3])- rrb'*f;
       f=f+mu*e*rrb;
       y(cnt)=rrb'*f;
       timer2=timer2+1;
       cnt=cnt+1;
    end
    timer2=0;
end

yflip=quantalph(y,[-3,-1,1,3])';
%% DEBUG - equalizer
if dm==1
    % EqualizerTest.m test routine 
    figure(11),plot(xs,'.');
    figure(12),plot(y,'.');
end
%% Correlator 2
framlent=784;
preamble=letters2pam('A0Oh well whatever Nevermind');
prelent=length(preamble);   
xsc=xcorr(preamble,y);                % do cross correlation
[C,ind]=max(abs(xsc));                % location of largest correlation...
headstart=length(y)-ind+1;            % ...gives place where header starts
headstart=mod(headstart,framlent)+framlent;

numpre=0;
g=length(y)-headstart;
while g>=0
    g=g-framlent;
    numpre=numpre+1;
end
%% DEBUG - Correlator 2
if dm==1
    subplot(3,1,1), stem(preamble)       % plot header
    title('Header')
    subplot(3,1,2), stem(y)             % plot data sequence
    title('Data with embedded header')
    subplot(3,1,3), stem(abs(xsc))            % plot correlation
    title('Correlation of header with data')
end
%% Decoding
savedheadstart=headstart;

newSlot = [];
oldSlot = [];
for count=1:numpre
    if ( headstart+prelent+1*datlent-1 <= length(yflip) )
        newSlot = yflip(headstart+prelent+0*datlent:headstart+prelent+1*datlent-1);
        message1 = [oldSlot newSlot];
        oldSlot = message1;
    else
        newSlot = yflip(headstart+prelent+0*datlent:length(yflip));
        message1 = [oldSlot newSlot];
        oldSlot = message1;
    end
    headstart=headstart+framlent;
end

headstart=savedheadstart;

newSlot = [];
oldSlot = [];
for count=1:numpre
    if ( headstart+prelent+2*datlent-1 <= length(yflip) )        
        newSlot = yflip(headstart+prelent+1*datlent:headstart+prelent+2*datlent-1);
        message2 = [oldSlot newSlot];
        oldSlot = message2;
    else
        newSlot = yflip(headstart+prelent+1*datlent:length(yflip));
        message2 = [oldSlot newSlot];
        oldSlot = message2;
    end
    headstart=headstart+framlent;
end

headstart=savedheadstart;

newSlot = [];
oldSlot = [];
for count=1:numpre
    if ( headstart+prelent+3*datlent-1 <= length(yflip) )
        newSlot = yflip(headstart+prelent+2*datlent:headstart+prelent+3*datlent-1);
        message3 = [oldSlot newSlot];
        oldSlot = message3;
    else
        newSlot = yflip(headstart+prelent+2*datlent:length(yflip));
        message3 = [oldSlot newSlot];
        oldSlot = message3;
    end
    headstart=headstart+framlent;
end

if desired_user==1
    decoded_text=pam2letters(message1)
end
if desired_user==2
    decoded_text=pam2letters(message2)
end
if desired_user==3
    decoded_text=pam2letters(message3)
end

end



